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A  Next- Generation  Model  of  the  Corona  and  Solar 

Wind:  Final  Report 


1  Introduction 

A  central  goal  of  solar  and  space  physics  is  to  understand  the  impacts  of  the  Sun  and  its 
activity  on  the  heliosphere,  particularly  the  space  environment  near  the  Earth.  While  coronal 
mass  ejections  (CMEs)  are  the  most  spectacular  example  of  this  influence,  the  structure  and 
dynamics  of  the  ambient  solar  corona  and  solar  wind  is  also  important.  Coronal  structure 
leads  to  the  partitioning  of  the  solar  wind  into  fast  and  slow  streams,  which  are  the  source  of 
recurrent  geomagnetic  activity.  The  geo-effectiveness  of  CMEs  is  in  part  determined  by  their 
interaction  with  the  ambient  wind.  The  connection  of  the  ambient  interplanetary  magnetic 
field  to  CME-related  shocks  and  impulsive  solar  flares  determines  where  solar  energetic  par¬ 
ticles  propagate.  If  we  are  to  make  substantial  progress  in  understanding  and  predicting  the 
Sun’s  space  weather  effects,  a  quantitative  description  of  the  the  solar  corona  and  solar  wind 
is  essential. 

The  goal  of  our  project,  which  has  been  supported  by  the  joint  AFOSR/NASA/NSF  part¬ 
nership  for  Collaborative  Space  Weather  Modeling,  is  to  improve  and  further  develop  a  time- 
dependent,  three-dimensional  MHD  model  of  the  solar  corona  and  solar  wind,  and  to  deliver 
versions  of  the  model  for  community  use.  Our  project  focuses  on  improving  the  underlying 
physics,  boundary  conditions,  and  usability  of  the  model  in  several  key  areas,  as  described 
in  this  report.  The  general  framework  for  our  coronal  and  heliospheric  solutions  is  described 
as  “CORHEL”  for  Corona-Heliosphere.  The  original  version  of  CORHEL,  developed  for  the 
CISM  (Center  for  Integrated  Space  Weather  Modeling)  project,  included  the  MAS  (Magneto¬ 
hydrodynamic  Algorithm  outside  a  Sphere)  code  for  coronal  computations  and  the  Enlil  code 
for  heliospheric  computations.  To  minimize  confusion  and  support  one  deliverable  product, 
we  have  retained  the  name  CORHEL  for  the  products  of  this  project.  As  a  result  of  our  joint 
support,  the  capabilities  of  CORHEL  have  been  greatly  improved  and  extended.  CORHEL 
now  delivers  solutions  to  the  community  via  four  customers  at  the  present  time:  the  CCMC, 
CISM,  our  own  modeling  web  site  (www.predsci.com),  and  AFRL  in  Albuquerque.  An  es¬ 
sential  aspect  of  our  modeling  approach  for  CORHEL  is  the  use  of  different  approximations 
for  the  coronal  solution,  and  to  allow  these  solutions  to  provide  the  inner  boundary  condi¬ 
tion  for  heliospheric  solutions.  A  flow  chart  of  the  key  components  of  CORHEL  is  given  in 
Appendix  A. 

The  great  utility  of  a  coronal  and  heliospheric  modeling  suite  is  the  ability  to  assist  in 
interpreting  both  remote  solar  and  in  situ  measurements  within  a  global  context.  To  illustrate 
this,  Figure  1  summarizes  the  state  of  the  solar  corona  and  inner  heliosphere  for  Carrington 
rotation  (CR)  1922,  which  occurred  from  April  24,  1997  through  May  21,  1997,  and  provided 
the  background  state  into  which  the  well-studied  May  12,  1997  CME  erupted  (Webb  et  al., 
2000;  Odstrcil  et  al.,  2004,  2005;  Linker  et  ah,  2007).  The  coronal  solution  on  the  left  was 
computed  using  our  most  sophisticated  thermodynamic  model  (see  section  2.1),  while  the 
heliospheric  solution  on  the  right  was  generated  using  a  simplified,  empirically  based  polytropic 
model.  In  the  following  sections  we  briefly  describe  the  capabilities  of  the  different  components 
of  CORHEL. 
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Structure  of  the  Corona  on  May  12, 1997 


Structure  of  the  Heliosphere  During  Carrington 
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Figure  1:  A  detailed  numerical  simulation  of  a  CME  can  capture  the  detailed  properties  of 
the  solar  source  active  region  (left),  including  measured  magnetic  fields  derived  from  magne¬ 
tograms,  as  well  as  the  structure  of  the  background  solar  wind  in  the  heliosphere  (right). 

2  The  Coronal  MHD  Model  MAS 

The  PSI  coronal  code  MAS  solves  the  MHD  equations  on  a  nonuniform  spherical  grid  that 
allows  us  to  concentrate  grid  points  in  regions  of  interest.  The  method  of  solution,  including 
the  boundary  conditions,  has  been  described  previously  (Mikic  &  Linker,  1994;  Linker  &  Mikic, 
1997;  Lionello  et  ah,  1999;  Mikic  et  al.,  1999;  Linker  et  al.,  2001;  Lionello  et  ah,  2001).  A 
description  of  the  equations  solved  in  the  MHD  model  is  given  in  Appendix  A.l.  To  compute 
a  coronal  solution,  the  MAS  model  integrates  these  time-dependent  equations  to  steady  state 
for  a  given  boundary  condition  defined  by  a  magnetic  map.  We  commonly  speak  of  two 
approximations  for  the  model.  The  polytropic  model  neglects  the  complicated  physics  of  the 
transition  region  by  setting  the  ratio  of  specific  heats  7  to  a  reduced  value.  The  model  retains 
an  energy  equation  (Eq.  (6)  in  Appendix  A.l  with  S  =  0)  but  we  refer  to  it  as  polytropic  to 
denote  this  limitation.  We  began  using  this  approximation  15  years  ago  in  our  first  eclipse 
calculation  (Mikic  &  Linker,  1996)  and  it  remains  useful;  it  is  still  used  in  many  contemporary 
MHD  models  (e.g.,  Cohen  et  al.,  2008;  Hayashi  et  al.,  2008).  A  more  advanced  treatment  of 
the  energy  equation  is  described  below. 

2.1  The  Thermodynamic  MHD  Model 

A  key  enhancement  in  the  CORHEL  model  is  the  inclusion  of  important  energy-transport  pro¬ 
cesses  (radiative  losses,  anisotropic  thermal  conduction,  and  coronal  heating)  in  the  transition 
region  and  solar  corona.  We  refer  to  this  as  the  “thermodynamic”  MHD  model.  This  more 
accurate  representation  of  energy  flow  in  the  MHD  model  allows  us  to  compute  simulated 
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Figure  2:  Comparison  between  the  emission  observed  on  August  27,  1996  in  SOHO/EIT  and 
Yohkoh/SXT  (Top)  and  a  coronal  MHD  simulation  with  energy  transport  terms  included  in 
the  energy  equation  (Bottom) .  The  images  from  the  calculation  were  synthesized  by  integrat¬ 
ing  the  appropriate  emission  kernels  along  the  line  of  sight.  Observations  and  simulations  are 
scaled  identically. 


EUV  and  X-ray  emission  as  is  observed  from  instruments  such  as  EIT  aboard  SOHO,  the 
SECCHI  EUV  imagers  aboard  the  STEREO  spacecraft,  the  Hinode  XRT  instrument,  and  the 
AIA  instrument  aboard  SDO.  Comparison  with  actual  emission  provides  powerful  constraints 
on  the  models. 

The  principal  uncertainty  in  the  thermodynamic  MHD  model  arises  from  the  unknown 
mechanism  of  coronal  heating.  Although  there  has  been  considerable  theoretical  activity  in 
understanding  coronal  heating  and  acceleration  of  the  solar  wind  (e.g.,  Piddington,  1956; 
Osterbrock,  1961;  Parker,  1972,  1994;  Roberts,  2000;  Velli,  2003;  Hollweg,  2003),  the  details 
are  not  understood. 

In  order  to  produce  models  that  are  consistent  with  observations,  we  have  developed  a 
parameterized  approach  for  describing  the  coronal  heating  (Lionello  et  ah,  2009).  Although 
this  limits  our  ability  to  address  theories  of  heating  in  any  detail,  it  does  provide  a  pragmatic 
compromise  of  being  able  to  include  the  chromosphere  and  transition  region  within  our  code 
and  restore  the  adiabatic  index,  7,  back  to  5/3.  Additionally,  comparisons  between  simulated 
emission  (EUV  and  soft  X-ray)  and  observations  provide  strong  constraints  on  the  free  pa¬ 
rameters  in  the  heating  model  (Lionello  et  al.,  2009).  We  are  also  developing  heating  models 
derived  from:  (1)  an  extension  of  robust  scaling  laws  (Rappazzo  et  ah,  2007)  based  on  Parker’s 
idea  of  shaking  and  tangling  of  magnetic  field  lines  (Parker,  1972)  to  include  shell  models  of 
turbulence  (Buchlin  &  Velli,  2007);  and  (2)  Wave/ turbulence- driven  (WTD)  models,  which 
attempt  to  provide  a  self-consistent  description  of  both  the  acceleration  and  heating  of  solar 
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Comparison  of  Simulated  and  Observed  Emission  for  Four  Different  Time  Periods 


Simulated  SOHO  EIT  and  Yohkon  SXT  Emission  [Log10DN/s] 


EIT171A  EIT195A  EIT284A  SXT  Al-Mg  Filter 


(c)  Observed  Emission  on  September  20, 2007  near  13:00UT  [Log10DN/s] 


Simulated  SOHO  EIT  and  Yohkon  SXT  Emission  [Log10DN/s] 


EIT171A  EIT195A  EIT284A  XRT  Al-Mesh  Filter 


Simulated  STEREO  A  EUVI  and  Hinode  XRT  Emission  [Log10DN/s] 
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Figure  3.  Comparison  of  simulated  and  observed  emission  for  different  time  periods  with  the  same  composite  heating  model:  (a)  May  11, 1997.  This  simulation 
was  used  as  the  background  for  simulations  of  the  May  12,  1997  CME.  (b)  May  13,  2005.  This  simulation  is  also  the  prelude  to  studies  of  the  May  2005  CME. 
Much  higher  field  strength  in  the  active  regions  yields  higher-than-observed  emission  with  this  heating  model,  (c)  September  20,  2007.  The  Sun  was  very  quiet 
during  this  time.  An  extended  coronal  hole  is  visible  in  the  observations  and  the  simulation,  (d)  July  19,  2008.  This  is  the  rotation  prior  to  the  August  8,  2008 
eclipse.  The  Sun  is  again  very  quiet.  The  three  left-most  images  are  from  EUVI  aboard  the  STEREO  A  spacecraft.  STEREO  A  views  a  different  portion  of  the 
Sun  than  the  Hinode  spacecraft,  where  the  XRT  image  on  the  right  is  taken  from. 


wind  plasma  through  the  combined  effects  of  wave  damping  and  turbulent  cascade  (Cranmer 
et  ah,  2007;  Cranmer,  2010).  These  approaches  will  have  fewer  free  parameters  and  may 
eventually  yield  more  accurate  solutions  than  a  purely  empirical  approach. 

In  Figure  2  we  compare  the  observed  images  with  those  synthesized  from  our  model  for 
August  27,  1996,  when  the  “Elephant’s  trunk”  equatorial  coronal  hole  was  visible.  To  make  a 
quantitative  comparison,  we  processed  SOHO/EIT  observations  using  the  SolarSoft  package 
to  produce  images  with  absolute  “DN/s”  [data  number/sec]  values,  in  which  the  images  are 
corrected  for  background  dark  current  subtraction,  degridding,  filter  normalization,  response 
sensitivity  correction,  and  exposure  time  normalization.  We  computed  the  emission  from  our 
MHD  solutions  by  integrating  the  emission  kernel,  f  n2ef(T)ds ,  over  the  line  of  sight  with  the 
appropriate  emission  and  instrument  response  function  f(T).  We  plot  the  log  of  the  emission 
with  identical  scales  for  the  observed  and  simulated  images.  Although  there  are  regions  in 
which  the  emission  is  dissimilar,  especially  on  small  length  scales  and  in  the  active  region,  the 
large  scale  features,  especially  the  coronal  hole  regions,  agree  reasonably  well. 

Figure  3  shows  comparisons  of  simulated  and  actual  emission  for  four  other  time  periods. 
Frame  (a)  shows  results  around  the  time  of  the  May  12,  1997  CME,  the  same  calculation 
illustrated  in  Figure  1.  Figure  3  illustrates  that  a  thermodynamic  model,  with  empirically- 
based  heating,  is  capable  of  reproducing  the  essential  features  of  the  observations  for  a  number 
of  time  periods.  Thus,  in  spite  of  the  fact  that  the  physical  processes  that  heat  the  corona 
remain  unknown,  we  have  captured,  at  least  to  a  first  approximation,  their  effects  within  our 
model. 


3  The  WSA  Model 

Support  from  this  program  has  allowed  us  to  includes  a  version  of  the  potential  field  source- 
surface/Wang-Sheeley-Arge  (WSA)  model  in  CORHEL  to  provide  quick-look  solutions.  The 
well-known  potential  field  source-surface  model  (PFSS)  is  a  potential  (current-free)  approx¬ 
imation  to  the  coronal  magnetic  field  obtained  by  solving  Laplace’s  equation  (V2<L  =  0) 
assuming  that  the  magnetic  field  becomes  radial  at  a  certain  height  (usually  2.5/?o,  where 
Ro  is  the  solar  radius).  A  variant  of  this  model  solves  for  the  field  above  the  source  surface, 
under  the  assumption  that  regions  of  opposite  magnetic  polarity  are  separated  by  a  current 
sheet  (potential  field  current  sheet  or  PFCS  model).  The  WSA  model  empirically  estimates 
the  solar  wind  speed  based  on  a  combination  of  the  expansion  factor  of  the  magnetic  field  and 
the  distance  of  magnetic  footpoints  from  coronal  hole  boundaries.  Zoran  Mikic  at  PSI  and 
Nick  Arge  at  AFRL  have  implemented  the  PFSS,  PFCS,  and  WSA  models  using  the  modern 
numerical  potential  solver  from  the  MAS  code.  Figure  4  shows  an  example  of  a  WSA  solution 
from  our  model. 

Incorporating  the  WSA  model  in  CORHEL  provides  a  number  of  advantages.  Solutions 
depend  sensitively  on  the  input  magnetic  maps.  Very  often,  comparisons  between  models 
really  turn  out  to  be  comparisons  between  how  the  input  data  was  processed.  CORHEL  can 
provide  WSA  and  MHD  solutions  using  the  same  magnetic  map  and  grid  so  that  different 
approaches  to  coronal  modeling  can  be  directly  compared  using  the  same  map.  Second,  as  the 
magnetic  maps  are  so  important,  the  WSA  solutions  allow  one  to  quickly  examine  the  results 
for  different  observatories  and  different  map  filtering.  Third,  this  new  version  of  WSA  allows 
the  development  of  higher  resolution  potential  field  solutions  than  is  possible  with  the 
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Figure  4.  (a)  Open  (black)  and  closed  (gray) 
magnetic  field  regions  for  the  potential  field 
source-surface  model.  This  model  was 
computed  using  the  numerical  potential  field 
solver  from  the  MAS  MHD  code,  (b)  Estimate 
of  the  solar  wind  speed  from  the  WSA  empiri¬ 
cal  formula  based  on  the  magnetic  field 
expansion  factor  and  distance  of  footpoints 
from  coronal  holes.  The  implementation  of 
this  capability  using  the  numerical  potential 
solver  allows  for  quick-look  solutions  prior  to 
running  a  full  MHD  model,  and  direct 
comparisons  between  empirical  MHD  and 
models. 
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Figure  5.  Heliospheric  solutions  computed  out  to  5  AU for  CR2068  (March  19-April  16,  2008,  the  WHI  inter¬ 
val).  Both  solutions  were  computed  with  the  MAS  heliospheric  model.  Velocity  color  maps  in  a  meridional 
slice  and  the  equatorial  plane  are  shown,  together  with  selected  magnetic  field  lines.  The  case  on  the  left  used 
the  WSA  model  for  the  interior  coronal  solution  and  specification  of  the  inner  heliospheric  boundary  condition. 
The  case  on  the  right  used  the  poly  tropic  MAS  coronal  model.  The  stream  interfaces  are  sharper  and  the  high¬ 
speed  streams  are  faster  when  the  MHD  model  is  used  for  the  coronal  solution. 
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standard  polynomial-based  potential  solvers,  as  well  as  nonuniform  resolution  (e.g.  high  res¬ 
olution  in  a  particular  active  region).  This  opens  the  possibility  that  further  experimentation 
could  improve  these  empirical  solutions. 


4  Heliospheric  Models 
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Figure  6:  Comparison  of  model  results  (red)  with  observations  (blue)  for  several  time  periods: 
(a)  Ulysses  rapid  latitude  scan;  and  (b)  Carrington  rotation  1913  (part  of  the  “Whole  Sun 
Month”  interval),  when  Ulysses  was  located  at  ~  28° N  heliographic  latitude,  at  a  distance  of 
~  4.3  AU  from  the  Sun,  and  on  the  opposite  side  of  the  Sun  from  Earth. 

The  three  forms  of  coronal  solutions  described  in  sections  2  and  3  can  be  coupled  to  helio¬ 
spheric  models,  with  the  results  at  the  outer  boundary  of  the  coronal  calculation  supplying  the 
inner  boundary  conditions  for  the  heliospheric  simulation.  We  briefly  describe  this  coupling 
in  Appendix  A. 2.  We  note  here  that  both  the  WSA  model  and  the  polytropic  MHD  model 
require  an  ad  hoc  prescription  for  the  boundary  velocity  to  achieve  realistic  solutions.  Even 
when  a  steady-state  coronal  solution  is  supplied  to  the  heliospheric  model,  the  heliospheric  cal¬ 
culation  is  performed  in  the  inertial  frame  and  so  becomes  a  time-dependent  calculation.  Two 
different  heliospheric  MHD  models  are  now  available  in  CORHEL;  the  Enlil  model  developed 
by  Dusan  Odstrcil,  and  a  heliospheric  version  of  the  MAS  model.  Figure  5  shows  an  example 
of  heliospheric  solutions  using  both  the  WSA  and  the  polytropic  MAS  model  for  the  coronal 
solution.  For  both  cases  the  heliospheric  MAS  model  was  used  for  the  outer  solution.  Note 
that  the  high-speed  streams  are  faster  and  the  stream  interfaces  more  sharply  pronounced 
when  the  MAS  model  is  used  for  the  coronal  solution.  This  appears  to  be  a  typical  result 
regardless  of  which  heliospheric  code  is  used. 

Our  team  at  PSI  have  used  global  MHD  models  to  interpret  the  observed  properties  of  the 
quasi-steady  corona  and  ambient  solar  wind  (Riley  et  al.,  1996,  2001a, b,  2002,  2003;  Riley, 
2007)  for  a  number  of  years.  In  general,  these  models  can  reproduce  the  essential  features  of 
the  observations.  In  Figure  6,  for  example,  we  compare  model  results  (computed  by  “flying” 
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the  spacecraft  trajectories  through  the  simulation  results)  with  Ulysses  observations  during: 
(a)  the  so-called  rapid  latitude  scan;  and  (b)  CR  1913,  at  which  time  Ulysses  was  returning 
to  lower  latitudes  and  was  located  at  ~  28° N  at  a  distance  of  ~  4.3  AU,  and  on  the  opposite 
side  of  the  Sun  from  Earth.  Comparisons  like  these  demonstrate  that  the  model  can  often 
reproduce  the  basic  features  of  the  solar  wind  under  quiescent  conditions.  However,  the 
model  does  not  always  perform  so  well,  particularly  when  the  coronal  magnetic  field  evolves 
significantly  from  one  rotation  to  the  next.  We  believe  that  future  models  that  use  evolving 
boundary  conditions  could  yield  significantly  more  accurate  solutions  than  have  been  possible 
in  the  past.  The  boundary  conditions  could  be  updated  in  time  using  a  combination  of  new 
magnetograms  as  they  become  available  and  phototospheric  flux  evolution  models  (Arge  et  al., 
2010). 

For  computing  heliospheric  solutions,  we  have  developed  two  complementary  approaches. 
In  the  simpler,  ad  hoc  technique  (illustrated  above),  we  use  the  structure  of  the  coronal 
magnetic  field  to  derive  the  radial  velocity  boundary  condition  for  the  inner  boundary  of  the 
heliospheric  model.  It  is  based  on  the  idea  that  the  slow  solar  wind  originates  at  the  boundary 
between  open  and  closed  field  lines.  We  use  the  computed  magnetic  field  from  the  coronal 
solution  directly,  and  infer  the  remaining  plasma  quantities  (density  and  temperature)  by 
assuming  momentum  flux  conservation  and  pressure  balance  over  the  sphere  defining  the  inner 
boundary.  The  second,  more  self-consistent  approach  is  to  run  the  heliospheric  model  directly 
using  all  of  the  magnetic  and  plasma  variables  computed  as  part  of  the  coronal  solution.  So  far, 
the  ad  hoc  solutions  tend  to  better  match  in  situ  data  at  Earth  and  Ulysses,  primarily  because 
the  transition  from  slow  to  fast  wind  is  too  gradual  with  our  present  heating/acceleration 
model.  Eventually,  we  expect  that  as  the  physics  contained  within  the  coronal  model  becomes 
more  sophisticated,  and  the  remaining  free  parameters  become  better  constrained,  the  quality 
of  the  self-consistently-derived  heliospheric  solutions  will  surpass  the  ad  hoc  results. 


Longitude  Longitude 


Figure  7:  Comparison  of  synoptic  maps  (as  a  function  of  longitude  and  longitude,  for  four 
solar  observatories.  Synoptic  Optical  Long-term  Investigations  of  the  Sun  (SOLIS),  Michel- 
son  Doppler  Imager  (MDI);  Global  Oscillation  Network  Group  (GONG);  and  Wilcox  Solar 
Observatory  (WSO).  The  Wang  correction  factor  was  applied  to  the  WSO  data. 
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5  Observationally-Derived  Boundary  Conditions 
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Figure  8:  Predicted  CORHEL  solar  wind  velocity  for  Carrington  Rotation  2060  (August  14  - 
September  10,  2007)  using  boundary  conditions  derived  from  (top  to  bottom):  SOLIS,  MDI 
aboard  SOHO,  NSO’s  GONG,  and  Wilcox  Solar  Observatory  (maps  shown  in  Fig.  7).  The 
frames  on  the  left  show  the  velocity  mapped  on  to  a  1  AU  sphere,  viewed  on  day  241  (August 
29,  2007)  from  Earth.  The  approximate  locations  of  the  STEREO  A,  B,  and  ACE  spacecraft 
are  shown.  The  position  of  the  heliospheric  current  sheet  for  the  different  solutions  is  marked 
as  a  gray  surface  in  these  frames.  The  frames  on  the  right  show  comparisons  of  the  predicted 
flow  speed  with  the  STEREO  A,  B,  and  ACE  spacecraft. 

The  Sun’s  magnetic  field  is  a  vital  ingredient  to  any  predictive  model  of  the  corona  and 
solar  wind.  Full-disk  measurements  of  the  line-of-sight  component  of  the  magnetic  field  in 
the  photosphere  have  been  available  for  almost  40  years.  As  no  consensus  has  emerged  as 
to  whether  a  specific  observatory  has  the  most  accurate  maps,  we  have  designed  CORHEL 
to  allow  the  user  to  choose  the  the  input  data  from  a  range  of  instruments.  The  user  can 
select  from  6  different  magnetographs  for  the  input  boundary  data:  Mount  Wilson  Observa¬ 
tory  (MWO);  the  National  Solar  Observatory  Vacuum  telescope  at  Kitt  Peak  (NSO/KP,  for 
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dates  prior  to  September,  2003);  the  NSO  vector  spectro-magnetograph  (NSO/SOLIS,  after 
September,  2003);  NSO/GONG,  the  Michelson  Doppler  Imager  aboard  the  SOHO  spacecraft 
(MDI),  and  Wilcox  Solar  Observatory  (WSO).  Magnetograms  are  now  available  from  the  He- 
lioseismic  and  and  Magnetic  Imager  (HMI)  aboard  the  SDO  spacecraft  launched  in  February, 
2010,  CORHEL  will  start  using  HMI  magnetic  maps  when  they  become  routinely  available. 

The  magnetic  maps  supplied  from  different  observatories  usually  agree  very  well  qual¬ 
itatively,  but  often  differ  quantitatively  These  differences  are  not  trivial  and  can  lead  to 
important  differences  in  the  solutions.  As  an  illustration,  Figure  7  summarizes  synoptic  maps 
as  a  function  of  longitude  and  latitude  from  four  solar  observatories  for  Carrington  rotation 
2060.  The  SOLIS,  MDI,  and  GONG  maps  were  interpolated  to  a  lower  resolution  180  x  90 
resolution;  the  WSO  is  displayed  at  nonuniform  36  x  73  resolution  (essentially  the  standard 
resolution  that  is  available).  The  structures  in  the  maps  are  similar,  but  there  are  quantitative 
differences  in  the  magnitude  of  the  fields.  The  global  structure  of  the  velocity  predicted  by 
the  MHD  models  using  these  different  maps  is  also  qualitatively  similar,  as  can  be  seen  in 
the  left  hand  frames  of  Figure  8.  However,  the  velocity  predicted  for  the  STEREO  and  ACE 
spacecraft  can  differ  significantly  when  different  maps  are  used,  as  shown  in  the  right  hand 
frames  of  Figure  8.  Because  Earth  sits  in  the  slow  wind  band  associated  with  the  coronal 
streamer  belt,  the  solar  wind  parameters  are  sensitive  to  quantitative  differences  in  the  maps, 
which  can  strongly  influence  the  accuracy  of  the  modeled  solutions  when  compared  with  1  AU 
measurements.  Much  of  this  difference  may  occur  because  of  differing  estimates  for  the  polar 
fields.  Our  results  indicate  that  the  1  AU  results  are  also  likely  to  be  sensitive  to  changes  in 
photospheric  field  that  occur  over  the  course  of  a  solar  rotation. 


6  Summary 

CORHEL  is  a  coupled  set  of  models  and  tools  for  quantitatively  modeling  the  ambient  solar 
corona  and  solar  wind  in  various  approximations.  Support  from  the  AFOSR/NASA/NSF 
partnership  for  Collaborative  Space  Weather  Modeling  has  allowed  us  to  greatly  enhance  the 
capabilities  of  CORHEL.  CORHEL  is  available  for  community  use  at  the  CCMC,  and  is  also 
being  run  by  CISM  and  AFRL.  We  are  presently  working  with  the  CCMC  to  enhance  the 
user  interface  for  selecting  model  parameters  and  diagnostic  outputs.  We  have  used  CORHEL 
and/or  its  component  models  in  a  number  of  studies,  as  shown  in  the  publications  section. 
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A  Appendix:  Model  Description 

Figure  9  shows  the  present  functionality  of  CORHEL.  In  sections  A.1-A.3  we  describe  some 
details  of  CORHEL,  including  the  equations,  boundary  and  initial  conditions  for  the  coronal 
model,  and  coupling  of  coronal  and  heliospheric  solutions. 


Magnetic  Maps:  MDI,  MWO,  NSO/KP,  NSO/GONG,  NSO/SOLIS,  WSO 
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Figure  9:  Flow  chart  depicting  the  present  state  of  CORHEL.  Green  boxes  indicate  the 
primary  input:  data  derived  from  measured  solar  magnetic  fields.  The  coronal  model  en¬ 
compasses  three  choices:  WSA  runs  for  immediate  answers  from  a  baseline  model,  polytropic 
MHD  for  quick  turnaround  runs,  and  thermodynamic  MHD  for  the  most  physically  realistic 
description.  Red  indicates  coronal  model  output,  including  simulated  emission  and  white  light 
images.  Coronal  models  feed  into  the  heliospheric  model  (two  choices)  using  either  empirical 
prescriptions  for  the  velocity,  density,  and  temperature  or  directly  driven  by  the  outputs  from 
the  thermodynamic  coronal  model. 


A.l  Coronal  Modeling:  Equations 

The  PSI  coronal  code  MAS  (Magnetohydrodynamic  Algorithm  outside  a  Sphere)  solves  the 
following  equations  in  spherical  coordinates: 


47 r  .  . 

V  x  B  =  — J, 

r 

(1) 

1  <9B 

V  x  E=  —  — , 
c  at 

(2) 

E  +  -v  x  B  =  rjJ, 

(3) 
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(4) 


!  +  V.(pv)=0, 

=  x  B  -  Vp  -  Vpw  +  pg  +  V  •  (z/pVv),  (5) 

T^r(f +v'VT)  =  -rv'v+s'  (6) 

5  =  2^T  (_V  '  q  "  n +  H  +  Hd  +  D )  (7) 

where  B  is  the  magnetic  field,  J  is  the  current  density,  E  is  the  electric  field,  p,  v,  p,  and 
T  are  the  plasma  mass  density,  velocity,  pressure,  and  temperature,  respectively,  ne  is  the 
electron  density,  and  pw  is  the  wave  pressure  and  represents  the  acceleration  due  to  waves. 
The  gravitational  acceleration  is  g  =  —g0rRl/r2,  R0  is  the  solar  radius,  7  =  5/3  is  the 
ratio  of  specific  heats,  p  is  the  resistivity,  and  u  is  the  viscosity.  The  plasma  pressure  is 
p  =  (ne  +  np)kT,  where  for  a  hydrogen  plasma,  ne  =  np ,  and  p  =  mpnp,  where  rnp  is  the 
proton  mass.  In  practice,  the  vector  potential  A  is  used  to  advance  the  equations.  These 
equations  are  solved  on  nonuniform  meshes  that  allow  us  to  concentrate  grid  points  in  regions 
of  interest.  The  method  of  solution,  including  the  boundary  conditions,  has  been  described 
previously  (Mikic  &  Linker,  1994;  Linker  &  Mikic,  1997;  Lionello  et  ah,  1999;  Mikic  et  al., 
1999;  Linker  et  al.,  2001;  Lionello  et  ah,  2001).  The  equations  are  solved  on  staggered  meshes, 
which  facilitates  the  implementation  of  boundary  conditions  and  enforces  the  vector  identifies 
V  ■  Vx  =  V  x  V  =  0  to  round-off  error.  This  approach  ensures  that  V  •  B  =  0  in  the 
calculation. 

In  the  energy  equation  (6)— (7) ,  H  is  the  coronal  heating  source,  H,i  =  ?/J2  +  z/\7v :  Vv  is  the 
heating  due  to  resistive  and  viscous  dissipation,  D  is  the  heating  due  to  dissipation  of  Alfven 
waves,  ne  and  np  are  the  electron  and  proton  density,  and  Q(T)  is  the  radiation  loss  function 
(e.g.,  Rosner  et  al.,  1978;  Athay,  1986).  In  the  collisional  regime  (below  ~  10.Ro),  the  heat 
flux  is  given  by  q  =  —  /cybb-VT,  where  b  is  the  unit  vector  along  B,  and  =  9  x  10-7T5//2  is 
the  Spitzer  value  of  the  parallel  thermal  conductivity  (in  cgs  units) .  In  the  collisionless  regime 
(beyond  ~  10.Ro)>  the  heat  flux  is  given  by  q  =  anekTv,  where  a  is  a  dimensionless  parameter 
of  order  1  (Hollweg,  1978).  The  (unknown)  coronal  heating  source  H  is  a  parameterized 
function,  as  discussed  in  section  2.1. 

Since  the  acceleration  of  the  solar  wind  by  Alfven  waves  occurs  on  spatial  and  time  scales 
below  those  of  our  global  numerical  model,  the  wave  pressure  pw  is  evolved  using  an  equation 
for  the  time-space  averaged  Alfven  wave  energy  density  e.  In  the  present  version  of  the  model, 
we  use  the  WKB  approach  (Jacques,  1977;  Usmanov  et  al.,  2000): 


<9v 

dt 


+  v-Vv 


—  +  V  •  F  =  v  •  Vpw  -  D,  (8) 

where  F  =  (|v  +  va)s  is  the  Alfven  wave  energy  flux,  va  =  B /\J An p  is  the  Alfven  speed,  and 
Pw  =  \e-  The  Alfven  wave  velocity  is  va  =  ±674 ;  we  transport  two  Alfven  wave  fields  (waves 
parallel  and  antiparallel  to  B),  which  are  combined  to  give  e.  The  Alfven  wave  energy  density 
e  is  related  to  the  space-time  average  of  the  fluctuating  component  of  the  magnetic  field  5B 
by  e  =  (SB2)  / 47 r.  The  dissipation  term  D  expresses  the  nonlinear  dissipation  of  Alfven  waves 
in  interplanetary  space  and  is  modeled  phenomenologically  (Hollweg,  1978).  We  are  presently 
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experimenting  with  a  more  fundamental  wave-driven  turbulence  model  (Cranmer  et  al.,  2007; 
Cranmer,  2010). 

The  simplified  polytropic  model  is  obtained  by  setting  pw  =  0  in  Eq.  (5),  S  =  0  in  Eq. 
(7),  and  7  close  to  1  (e.g.,  7  =  1.05).  The  polytropic  model  is  popular  because  it  reduces  the 
numerical  requirements  for  3D  coronal  calculations. 

The  MHD  equations  in  the  corona  are  very  stiff,  since  they  span  multiple  time  and  length 
scales.  In  a  typical  coronal  calculation,  one  can  have  mesh  cells  as  small  as  ~  100  km,  and  as 
large  as  a  few  solar  radii  (~  106 km).  The  Alfven  and  sound  speeds  can  range  from  ~20km/s 
to  ^10,000 km/s.  The  Lundquist  number  can  be  ^105.  For  these  equations,  implicit  methods 
are  crucial.  Typically,  we  exceed  the  wave  Courant  number  (which  would  limit  the  time  step 
in  an  explicit  calculation)  by  factors  of  50-100,  and  we  greatly  exceed  the  time  steps  needed 
for  stable  advancement  of  diffusive  terms  (thermal  conductivity,  resistivity,  and  viscosity) .  We 
use  a  semi-implicit  technique  to  provide  unconditional  stability  to  Alfven  and  sound  waves 
(Mikic  et  ah,  1988;  Harned  &  Kerner,  1985;  Harned  &  Schnack,  1986;  Schnack  et  ah,  1987), 
and  fully  implicit  advancement  of  diffusive  terms.  The  very  large  sparse  matrices  arising  from 
these  methods  are  inverted  using  an  iterative  preconditioned  Conjugate  Gradient  method. 

The  MAS  code  is  written  in  FORTRAN95  using  the  Message  Passing  Interface  (MPI)  with 
three-dimensional  domain  decomposition.  The  code  is  portable  and  runs  on  many  systems 
(including  Mac  OS  X  and  various  flavors  of  Linux  that  are  used  on  NASA’s  and  NSF’s  mas¬ 
sively  parallel  computer  architectures).  We  are  able  to  perform  runs  with  tens  of  millions  of 
mesh  points  on  hundreds  to  thousands  of  processors. 

A. 2  Coronal  Modeling:  Initial  Conditions,  Boundary  Conditions, 
and  Solution  Procedure 

The  initial  and  boundary  conditions  for  our  models  of  the  solar  corona  are  described  in  detail 
by  Linker  &  Mikic  (1997).  The  key  boundary  condition  required  from  observations  is  the 
radial  magnetic  field  at  the  solar  surface,  Br0.  This  must  be  supplemented  with  conditions  on 
the  plasma  temperature  and  density  at  r  =  R0.  In  polytropic  models,  these  are  coronal  values 
typically  in  the  range  of  To  =  1.4  —  2  x  106  K  and  no  =  1  —  4  x  108  cm-3.  The  choice  of  values 
influences  the  relative  size  of  the  closed  and  open  field  regions  and  the  speed  of  the  solar  wind. 
In  thermodynamic  models,  T0  =  20,000  K  and  n0  =  2  x  1012cm-3,  representative  of  the  upper 
chromosphere.  In  these  models,  the  choices  for  coronal  heating  determine  the  properties  of 
the  solutions,  and  the  exact  choice  of  T0  and  n0  are  not  crucial  as  long  as  n0  is  large  enough 
to  maintain  a  chromosphere  in  the  presence  of  the  specified  heating.  For  both  polytropic 
and  thermodynamic  models,  the  velocity  parallel  to  the  magnetic  field  cannot  be  specified 
in  a  well-posed  problem;  it  is  determined  self-consistently  at  r  =  R0  from  the  characteristic 
equations  along  B.  At  the  upper  boundary,  where  the  flow  is  supersonic  and  super- Alfvenic, 
all  quantities  are  computed  from  characteristic  equations. 

For  the  initial  condition,  we  start  by  computing  a  potential  magnetic  field  in  the  corona 
that  matches  Br0  at  the  solar  surface.  We  impose  a  spherically  symmetric  solar  wind  solution 
and  integrate  the  time-dependent  MHD  equations  in  time  until  the  solution  settles  down  to 
an  equilibrium.  Helmet  streamers  with  closed  field  lines  form,  surrounded  by  open  field  lines 
along  which  the  solar  wind  flows  outward. 
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A. 3  Heliospheric  Solutions 

To  compute  heliospheric  solutions,  the  results  at  the  outer  boundary  of  the  coronal  calculation 
must  supply  the  inner  boundary  conditions  for  the  heliospheric  simulation.  The  interface 
between  the  coronal  and  heliospheric  regions  is  located  in  the  super-critical  flow  region,  usually 
between  20.Ro  and  30Rq .  This  makes  the  coupling  procedure  relatively  straightforward,  since 
information  is  passed  only  one  way:  from  the  coronal  model  to  the  heliospheric  model.  Because 
the  physics  of  the  coupling  is  not  complicated,  a  sophisticated  coupling  architecture  is  not 
required.  The  coupling  is  accomplished  by  reading  and  writing  files  in  a  specified  format.  The 
procedure  is  very  modular  because  any  coronal  or  heliospheric  model  can  be  substituted;  the 
models  just  have  to  be  capable  of  reading/ writing  files  in  the  designated  format. 

To  produce  realistic  solar  wind  solutions  in  the  inner  heliosphere,  the  WSA  model  uses  an 
empirical  prescription  to  set  the  plasma  properties  at  the  inner  boundary  of  the  heliospheric 
model.  Polytropic  coronal  models  lack  the  necessary  physics  to  predict  realistic  contrast 
between  the  fast  and  slow  wind.  Therefore,  we  also  use  an  empirical  technique  (Riley  et  al., 
2001a)  to  specify  the  inner  boundary  conditions  for  the  velocity;  this  specification  has  some 
similarities  to  the  WSA  specification.  An  advantage  of  the  thermodynamic  MHD  model 
(section  2.1)  is  that  such  a  specification  is  not  necessary.  Given  a  specified  model  for  the 
coronal  heating  and  solar  wind  acceleration,  fast  and  slow  wind  regions  are  produced  self- 
consistently  in  the  coronal  model  and  can  be  fed  directly  to  the  heliospheric  model.  Even  in 
the  case  of  a  steady-state  coronal  solution,  the  heliospheric  solutions  are  time-dependent.  This 
is  because  the  steady-state  coronal  input  co-rotates  in  the  inertial  frame  of  the  heliospheric 
calculation. 

Both  the  MAS  and  Enlil  heliospheric  models  in  CORHEL  solve  ideal  or  nearly  ideal  ver¬ 
sions  of  Eqs.  (1-6).  The  outer  boundary  conditions  are  again  in  the  regime  of  supersonic  and 
super- Alfvenic  flow  and  are  computed  using  characteristic  equations. 
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